%%  粒子群
clc,clear;
close;
warning off;
% 参数
numParticles=50; %最大粒子数
maxIter=100;     %迭代次数
w = 0.5;          % 惯性权重
c1 = 1;         % 个体学习因子
c2 = 1;         % 全局学习因子
alapa=2;      %莱维飞行幂律分布的指数(0~2)

x=zeros(1,3);
numVars = length(x);    % 变量的数量
[bestx,bestf,bd]=lfpso(maxIter,w,c1,c2,numParticles,numVars,alapa);
disp("Best Solvution");
disp(bestx);
disp("Best Function");
disp(bestf);
tu(bd,maxIter);

%% 作图
function tu(bd,maxiter)
    plot(1:maxiter,bd,'r');
    xlabel('迭代次数')
    ylabel('适应值')
    title('过程')
end


%% PSO
function [bestx,bestf,bd]=lfpso(maxIter,w,c1,c2,numParticles,numVars,alpha0)
% 初始化粒子位置和速度
particlePositions = rand(numParticles, numVars) .* 10; % 初始化在0到10之间
particleVelocities = zeros(numParticles, numVars); % 初始速度为零
pbestx = particlePositions; % 每个粒子的历史最优位置
pbestf = inf(numParticles, 1); % 每个粒子的历史最优值
% 初始化全局最优
[bestf, bestParticleIndex] = min(pbestf);
bestx = pbestx(bestParticleIndex, :);
% PSO主循环
for iter = 1:maxIter
    for i = 1:numParticles
        % 更新粒子速度和位置
        for i = 1:numParticles
            r1 = rand(1, numVars);
            r2 = rand(1, numVars);
            particleVelocities(i, :) =w.*particleVelocities(i, :)+c1*r1.*((pbestx(i,:)+bestx)./2-particlePositions(i, :))+c2*r2.*((pbestx(i,:)-bestx)./2-particlePositions(i, :));
            particlePositions(i, :) = particlePositions(i, :) + particleVelocities(i, :);
        end
        % 检查粒子位置的可行性
        if yuesu3(particlePositions(i,:))
            currentValue =fun3(particlePositions(i,:)); % 计算目标函数值
        else
            currentValue = inf; % 不可行解
        end
        
        % 更新个人最优
        if currentValue <pbestf(i)
            pbestf(i) = currentValue;
            pbestx(i, :) = particlePositions(i, :);
        end
        % 更新全局最优
         [minValue, bestParticleIndex] = min(pbestf);
        if minValue < bestf
            bestf = minValue;
            bestx = pbestx(bestParticleIndex, :);
        end
    end
    % 应用莱维飞行
    for i = 1:numParticles
      step_length = levy_flight(alpha0); % 使用alpha = 1.5
      particlePositions(i, :) = particlePositions(i, :) + step_length * randn(1, numVars);
        % 检查粒子位置的可行性
        if yuesu3(particlePositions(i,:))
            currentValue1 =fun3(particlePositions(i,:)); % 计算目标函数值
        else
            currentValue1 = inf; % 不可行解
        end
        % 更新个人最优
        if currentValue1 <pbestf(i)
            pbestf(i) = currentValue1;
            pbestx(i, :) = particlePositions(i, :);
        end
        % 更新全局最优
         [minValue, bestParticleIndex] = min(pbestf);
        if minValue < bestf
            bestf = minValue;
            bestx = pbestx(bestParticleIndex, :);
        end
    end
     bd(iter)=bestf;
end
end
function h=yuesu3(x)
 h= ((x(1)>1) && (x(2) > 0)&&(x(1)+x(2)<10)&&(x(3)+x(2)>3));
end
function y=fun3(x)
%目标函数
    y=x(1)^2+2*x(1)*x(2)+3*x(2)^3+x(1)*x(2)^2-x(1)+x(3)^3;
end
% 莱维飞行生成函数
function L = levy_flight(beta)
    % 参数
    sigma_u = (gamma(1 + beta) * sin(pi * beta / 2) / (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2)))^(1 / beta);
    sigma_v = 1;

    % 生成莱维分布步长
    u = randn * sigma_u;
    v = randn * sigma_v;
    L = u / abs(v)^(1 / beta);
end
